/*
    Copyright 2010-2011 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/// \file
/// Classes for camera projections.

// $Id$

#ifndef __projection__
#define __projection__

#include "mcrx-types.h"


#include "blitz/numinquire.h"
#include "mcrx-units.h"
#include "constants.h"
#include "vecops.h"
#include "boost/shared_ptr.hpp"
#include <string>

namespace mcrx {
  class projection;
  class rectilinear_projection;
  class aitoff_projection;
}


namespace CCfits {
  class ExtHDU;
}


/** Abstract base class for camera projections. These map an incoming
    direction vector in a frame where the z-axis is normal to the
    image plane to normalized image positions (0,1) in x and y. */
class mcrx::projection {
public:

  /// Returns a string describing the projection type.
  virtual const std::string& type () const=0;
  /// Returns the total solid angle of the image.
  virtual T_float solid_angle () const=0;
  /** Returns the normalized image position (mapped to (0,1) resulting
      from projecting the position onto the image plane. If the
      position is not mapped onto the image, (NaN,NaN) is
      returned. Note that the mapping is always (0,1) even for
      non-unity aspect ratios. */
  virtual blitz::TinyVector<T_float, 2> project (vec3d pos) const=0;
  virtual void write_parameters(CCfits::ExtHDU&, const T_unit_map&) const=0;
  /** Returns the image scale in radians per normalized image position
      coordinate in the y-direction. (For small angle approximations,
      this is equal to the total field of view in radians.) Note that
      the NORMALIZED image scale in the x-direction depends on the
      projection, but the final image scale, once the differing
      numbers of pixels are taken into account, is always
      isotropic. */
  virtual T_float scale() const=0;
  /// Returns the aspect ratio of the image (pixels_x/pixels_y).
  virtual T_float aspect_ratio() const=0;

  static boost::shared_ptr<projection> 
  projection_factory(const std::string&, T_float, T_float);

  static boost::shared_ptr<projection>
  read_projection(CCfits::ExtHDU& hdu);

  /** "Deprojects" a normalized image position (0,1) to an outward
      direction. This is essentially the inverse of project(). */
  virtual vec3d deproject (const blitz::TinyVector<T_float, 2>&) const=0;
  /** Calculates dOmega/dA, the scale factor between solid angle of
      directions and the area on the normalized image position
      (0,1). This is not necessarily uniform unless a true
      area-preserving projection is being used. The input is the
      position on the normalized image plane. */
  virtual T_float area_scale(const blitz::TinyVector<T_float, 2>&) const=0;
};

/** Class for normal rectilinear projection cameras. This used to be
    the only alternative.  */
class mcrx::rectilinear_projection : 
  public mcrx::projection {
private:
  const T_float fov_; ///< Full field of view in radians in the y-dir.
  const T_float aspect_; ///< Aspect ratio (fov_x/fov_y).
  const vec3d lpn_; ///< Projection plane normal vector.
  const vec3d lpp_; ///< Projection plane center point.
  const vec3d ny_;  ///< y base vector in projection plane.
  const vec3d nx_;  ///< x base vector in projection plane.

  // dir_=zhat, so lpn=-zhat, nx_=xhat, ny_=yhat, right?

  friend boost::shared_ptr<projection>
  projection::read_projection(CCfits::ExtHDU& hdu);

  rectilinear_projection(CCfits::ExtHDU& hdu);

public:
  /// The type string identifying the projection
  static const std::string typestring_;

  /** Constructor which takes position and viewing direction as
    3-vectors. The "up" (y) direction is specified as another vector
    (which must not be parallel with the viewing direction.) In this
    case, the field of view is \b angular and in radians. */
  rectilinear_projection(T_float fov, T_float aspect) :
    fov_ (fov), aspect_(aspect),
    lpn_ (vec3d(0,0,-1)),
    //lpp_ (pos - lpn_/tan(fov_/2)),
    lpp_ (vec3d(0,0,1./tan(0.5*fov_))),

    // ny is the normalized projection of up orthogonal to lpn, standard
    // Gram-Schmidt stuff. (It is the responsibility of the user to make
    // sure this is nonsingular, but it does not have to be normalized)
    //ny_ ((up - dot(lpn_, up)*lpn_)/sqrt(dot((up - dot(lpn_, up)*lpn_),
    //(up - dot(lpn_, up)*lpn_)))),
    ny_(vec3d(0,1,0)),
    nx_ (cross (ny_, vec3d(0,0,-1))) 
  {
    assert (std::abs (dot (lpn_, lpn_) - 1) <1e-10);
    assert (std::abs (dot (nx_, nx_) - 1) <1e-10);
    assert (std::abs (dot (ny_, ny_) - 1) <1e-10);
    assert (std::abs (dot (lpn_, ny_) < 1e-10));
  };

  virtual const std::string& type () const { return typestring_; };
  virtual T_float solid_angle () const { return fov_*fov_; } ;
  virtual blitz::TinyVector<T_float, 2> project (vec3d pos) const;
  virtual void write_parameters(CCfits::ExtHDU&, const T_unit_map&) const;
  virtual T_float scale() const {return 2*tan(0.5*fov_); };
  virtual T_float aspect_ratio() const { return aspect_; };
  virtual vec3d deproject (const blitz::TinyVector<T_float, 2>&) const;
  virtual T_float area_scale(const blitz::TinyVector<T_float, 2>&) const;
};


/** Class for all-sky Hammer-Aitoff cameras. */
class mcrx::aitoff_projection : 
  public mcrx::projection {
private:
  // no data members necessary

public:
  /// The type string identifying the projection
  static const std::string typestring_;

  virtual const std::string& type () const { return typestring_;};
  virtual T_float solid_angle () const { return 4*constants::pi;} ;
  virtual blitz::TinyVector<T_float, 2> project (vec3d pos) const;
  virtual void write_parameters(CCfits::ExtHDU&, const T_unit_map&) const;
  /** Returns the image scale in radians per normalized image position
      coordinate in the y-direction. Currently not functional for
      Aitoff projection. \todo When we merge the camera response
      functionality from MLT this should be fixed. */
  virtual T_float scale() const {return 0;};
  /// Aspect ratio for a Hammer-Aitoff projection is fixed at 2.
  virtual T_float aspect_ratio() const { return 2; };
  virtual vec3d deproject (const blitz::TinyVector<T_float, 2>&) const;
  virtual T_float area_scale(const blitz::TinyVector<T_float, 2>&) const;
};

 
#endif
